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We present a code to construct initial data for binary neutron star systems in which the stars 
are rotating. Our code, based on a formalism developed by Tichy, allows for arbitrary rotation 
axes of the neutron stars and is able to achieve rotation rates near rotational breakup. We com¬ 
pute the neutron star angular momentum through quasi-local angular momentum integrals. When 
constructing irrotational binary neutron stars, we find a very small residual dimensionless spin of 
~ 2 X 10“'*. Evolutions of rotating neutron star binaries show that the magnitude of the stars’ an¬ 
gular momentum is conserved, and that the spin- and orbit-precession of the stars is well described 
by post-Newtonian approximation. We demonstrate that orbital eccentricity of the binary neutron 
stars can be controlled to ~ 0.1%. The neutron stars show quasi-normal mode oscillations at an 
amplitude which increases with the rotation rate of the stars. 

PACS numbers: 04.20EX, 04.25.dk, 04.30.Db, 04.40.Dg, 04.25.NX, 95.30sf 


I. INTRODUCTION 

Several known binary neutron star (BNS) systems 
will merge within a Hubble time due to inspiral driven 
by gravitational radiation [T], most notably the Hulse- 
Taylor pulsar [5]. Therefore, binary neutron stars consti¬ 
tute one of the prime targets for upcoming gravitational 
wave detectors like Advanced LIGO and Advanced 
Virgo mi- The neutron stars in known binary pul¬ 
sars have fairly long rotation periods [T]. The system 
J0737-3039 [7] contains the fastest known spinning neu¬ 
tron star in a binary with a rotation period of 22.7ms. 
This system will merge within ~ 10® years through grav¬ 
itational wave driven inspiral. Globular clusters con¬ 
tain a significant fraction of all known milli-second pul¬ 
sars [T], which through dynamic interactions, may form 
binaries 0 M- Gravitational wave driven inspiral re¬ 
duces nniiii] the initially high eccentricity of dynamical 
capture binaries. Given the presence of milli-second pul¬ 
sars in globular clusters, dynamically formed BNS may 
contain very rapidly spinning neutron stars with essen¬ 
tially arbitrary spin orientations. Presence of spin in BNS 
systems does influence the evolution of the binary. For 
instance, in order to avoid a loss in sensitivity in GW 
searches, one needs to account for the NS spin m- Fur¬ 
thermore, early BNS simulations m of irrotational and 
corotational BNS systems found that the spin of coro¬ 
tating BNS noticeably increased the size of accretion 
discs occurring during the merger of the two NS. The 
properties of accretion discs and unbound ejecta are inti¬ 
mately linked to electromagnetic and neutrino emission 
from merging compact object binaries |14j . Understand¬ 


ing the behavior of rotating BNS systems is therefore 
important to quantify the expected observational signa¬ 
tures from such systems. These considerations motivated 
a recent interest in the numerical modeling of rotating bi¬ 
nary neutron star systems during their last orbits and 
coalescence. Baumgarte and Shapiro [T^, Tichy [TH] . 
and East et al m presented formalisms for construct¬ 
ing BNS initial data for spinning neutron stars. Tichy 
proceeded to construct rotating BNS initial data [15] : 
and Ref. [19] studies short inspirals and mergers of BNS 
with rotation rates consistent with known binary neutron 
stars (i.e. a dimensionless angular momentum of each 
star X = S/M"^ < 0.05), and rotation axes aligned with 
the orbital angular momentum. Very recently, Dietrich 
et al. [20] presented a comprehensive study of BNS, in¬ 
cluding a simulation of a processing, merging BNS. East 
et. al. [H] investigate interactions of rotating neutron 
stars on highly eccentric orbit. Kastaun et al. [22] de¬ 
termine the maximum spin of the black hole remnant 
formed by the merger of two aligned spin rotating neu¬ 
tron stars. Tsatsin and Marronetti [S^ present initial 
data and evolutions for non-spinning, spin-aligned and 
anti-aligned data sets. 

Previous studies differ in the type of initial data used: 
Refs, [ini [m mi m] construct and utilize constraint- 
satisfying initial data, which also incorporates quasi¬ 
equilibrium of the binary system. Refs. [niEi] construct 
constraint-satisfying data based on individual TOV stars, 
without regard of preserving quasi-equilibrium in the re¬ 
sulting binary, but providing greatly enhanced flexibility 
in the type of configurations that can be studied, e.g. 
hyperbolic encounters. Refs. [5^ |13], finally, only ap- 


2 


proximately satisfy the constraint equations. Previous 
studies also differ in the rigor with which the neutron 
star angular momentum is measured. Ref. m merely 
discusses the neutron stars based on a rotational velocity 
oj® entering the initial data formalism (cf. our Eq. (48) 
below), whereas Refs. [T9l[2Tl[22] estimate the initial neu¬ 
tron star spin either based on single star models or based 
on the differences in binary neutron star initial data sets 
with and without rotation, and thus neglecting the im¬ 
pact of interactions in the binary. All these studies mea¬ 
sure the neutron star angular momentum in the initial 
data. Changes in the neutron star angular momentum 
that could happen during initial relaxation of the binary 
or during the subsequent evolution of the binary are not 
monitored. 

In this paper we study the construction of rotating bi¬ 
nary neutron star initial data and the evolution through 
the inspiral phase. We implement the constant rotational 
velocity (CRV) formalism developed by Tichy [18], and 
construct constraint satisfying BNS initial data sets with 
a wide variety of spin rates, as well as different spin di¬ 
rections. We apply quasi-local angular momentum tech¬ 
niques developed for black holes to our BNS initial data 
sets; the quasi-local spin indicates that we are able to 
construct BNS with dimensionless angular momentum 
exceeding 0.4. Evolving some of the constructed initial 
data sets through the inspiral phase, we demonstrate that 
we can control and reduce orbital eccentricity by an it¬ 
erative adjustment of initial data parameters controlling 
orbital frequency and radial velocity of the stars, both for 
non-precessing (i.e. aligned-spin binaries) and precessing 
binaries. When monitoring the quasi-angular momentum 
of the neutron stars during the inspiral, we find that its 
magnitude is conserved, and the spin-direction processes 
in a manner consistent with post-Newtonian predictions. 

This paper is organized as follows. Section [IT| describes 
the initial data formalism and our numerical code to solve 
for rotating BNS initial data. In Sec. we use this code 
to study a range of initial configurations, with a special 
emphasis on the behavior of the quasi-local spin diag¬ 
nostic. We evolve rotating BNS in Sec. EYl including 
a discussion of eccentricity removal, the behavior of the 
quasi-local spin diagnostics, and a comparison of the pre¬ 
cession dynamics to post-Newtonian theory. A discussion 
concludes the paper in Sec. |Vj In this paper, we work in 
units where G = c = Mq = 1. 


II. METHODOLOGY 
A. Formalism for irrotational binaries 

To start, we will review a formalism commonly used for 
the construction of initial data for system of irrotational 
binary neutron stars. We will then discuss how to build 
upon this formalism to construct initial data for neutron 
stars with arbitrary spins. 

We begin with the 3-1-1 decomposition of the space¬ 


time metric (see [24] for a review), 

ds^ = —a^dt^ jij (dx® -I- /3®dt) (da;-^ -I- P^dtj . (1) 

Here, a is the lapse function, /3® is the shift vector 
and is the 3-metric induced on a hypersurface S(t) of 
constant coordinate time t. In this decomposition, the 
unit normal vector n'® to E(t) and the tangent vector 
to the coordinate line t are related by 

= an^ + /3^, (2) 

with = (—a, 0,0,0) and (3^ = (0,/3®). The extrinsic 
curvature of S(t) is the symmetric tensor defined as 

fiu — ^u'^^jL ^1/7 ^^A(li^Oi) — (^) 

where is the extension of the 3-metric 

7 y- to the 4-dimensional spacetime, and is the 4- 
metric of that spacetime. By construction, = 0 

and we can restrict to the 3-dimensional tensor AT®-! 
defined on E x S. The extrinsic curvature AT®-’ is then 
divided into its trace K and trace-free part A®-^ : 

= A®^ + -7®^'a:. (4) 

3 

We treat the matter as a perfect fluid with stress- 
energy tensor 


= ip-\-P)Ufj,u^Pgij_^, (5) 

where p = po{l e) is the energy density, po the baryon 
density, e the specific internal energy, P the pressure, and 
Ufj_ the fluid’s 4-velocity. For the initial value problem, it 
is often convenient to consider the following projections 
of the stress tensor: 


E = 

(6) 


(7) 


(8) 

We then further decompose the metric according to 
the conformal transformation 

II 

A 

(9) 

Other quantities have the following conformal transfor¬ 
mations: 

E = 

(10) 

s = ^-^s, 

(11) 

J® = 

(12) 

II 

1 

o 

(13) 

a = 'I'®d. 

(14) 


A®-! is related to the shift and the time derivative of 
the conformal metric, Uij = dt^ij by 


A®^' = — 
2d 




( 15 ) 
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where L is the conformal longitudinal operator whose 
action on a vector 1 /® is 

(16) 

and V is the covariant derivative defined with respect to 
the conformal 3-metric 7 ^. 

In the 3+1 formalism, the Einstein equations are de¬ 
composed into a set of evolution equations for the metric 
variables as a function of t, and a set of constraint equa¬ 
tions on each hypersurface S(t). The initial data prob¬ 
lem consists in providing quantities g^Lv{to) and 
which satisfy the constraints on E(to) and represent ini¬ 
tial conditions with the desired physical properties (e.g. 
masses and spins of the objects, initial orbital frequency, 
eccentricity, etc.).We solve the constraint equations using 
the Extended Conformal Thin Sandwich (XCTS) formal¬ 
ism [ 25 ], in which the constraints take the form of five 
nonlinear coupled elliptic equations. The XCTS equa¬ 
tions can be written as 


2d 







= 0 , 


(17) 


V2^- 

8 12 

+ + 2Tr'S-^E = 0, (18) 
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ii?+^vI/4K2 + !vI;-84.ib 


+27r«'-2 (E + 2S) = {dtK - p^dkK) . (19) 


We solve these equations for the conformal factor 4', 
the densitized lapse d^'^ and the shift /3b E^ S and J® 
determine the matter content of the slice. The variables 
liE ~ K and dtK are freely chosen. 

If we work in a coordinate system corotating with 
the binary, utj = 0 and dtK = 0 are natural choices 
for a quasi-equilibrium configuration. Following earlier 
work [2H1428) . we also choose to use maximal slicing, 
K = 0, and a conformally flat metric, ■jij = 6ij. Max¬ 
imal slicing is a gauge choice that determines the loca¬ 
tion of the initial data hypersurface in the embedding 
space time. Conformal flatness is used for computational 
convenience; rotating black holes are known to be not 
conformally flat |29j . and so this simplifying assumption 
should be revisited in the future. 

In addition to solving these equations for the met¬ 
ric variables, we must impose some restrictions on the 
matter. In particular, the stars should be in a state 
of approximate hydrostatic equilibrium in the comoving 


frame. This involves solving the Euler equation and the 
continuity equation. For an irrotational binary, the first 
integral of the Euler equation leads to the condition 

ha^=C, ( 20 ) 

70 

where C is a constant, hereafter referred to as the Euler 
constant, the enthalpy h is defined as 


p 

h = 1 + e H -, 

Po 

(21) 

and we have introduced 


7 = 7n7o (1 - , 

(22) 

70 = (1 - 7ijt^ot^o) ' > 

(23) 

7„ = (1 - , 

(24) 

a- 

(25) 

The 3-velocity W is defined by 


= 7„(n'^ + C/^), 

(26) 

= 0. 

(27) 


The choice of t/*, which is unconstrained in this formal¬ 
ism, is an important component is determining the initial 
conditions in the neutron star. For irrotational binaries 
(non-spinning neutron stars), there exists a potential (j) 
such that 

W = ——^d,(j). (28) 

The continuity equation can then be written as a second- 
order elliptic equation for (j)' 

^V'^V^</)+(V^<(.)V^^ =0. (29) 


Under the assumption of the existence of an approx¬ 
imate helicoidal Killing vector ^ [301EU, this equation 
becomes 


Po 


- ^'^^didjcj) -5i7„ + 


a 


h 


dkfj) 


= 7 ■'diCpdjpo - dipo- 


(30) 


Another simple choice for t/* is to enforce corotation of 
the star, i.e. [/* = I/q. This would be the case if neutron 
star binaries were tidally locked. However, viscous forces 
in neutron stars are expected to be insufficient to impose 
tidal locking [32], and the neutron star spins probably 
remain close to their value at large orbital separations. 

Once we have obtained h from the metric and [/*, the 
other hydrodynamical variables can be recovered if we 
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close the system by the choice of an equation of state for 
cold neutron star matter in /3-equilibrium, P = P{po) 
and e = e(/ 0 o)- Throughout this work, we use a polytropic 
equation of state, P = kPq, with T = 2. The internal 
energy, epo, satisfies 


ePo = (31) 

The boundary conditions of our system of equations 
are quite simple. At the outer boundary of the com¬ 
putational domain (which we approximate as “infinity” 
and is in practice IO^^Mq), we require the metric to be 
Minkowski in the inertial frame, and so in the corotating 
frame we have 


f3 = flo X r + dor, 
a = 1, 

4- = l, 


(32) 

(33) 

(34) 


with Qq the initial orbital frequency of the binary and do 
the initial inspiral rate of the binary. We choose fio = 
( 0 , Ojl^o)) with rig and do as freely specifiable variables 
that determine the initial eccentricity of the binary. 

At the surface of each star, the boundary condition can 


be easily inferred from the po = 0 limit of equation (30): 


7 ■'dicpdjpo = - Oipo- 


(35) 


Finally, we discuss how a first guess for the orbital 
angular velocity rig can be obtained for a non-spinning 
system. The force balance equation at the centre of the 
NS is 


Vln/i = 0. (36) 

Neglecting any infall velocity, this condition guarantees 
that the binary is in a circular orbit. This is only an 
approximation as there is really some infall velocity, but 
this still leads to low eccentricity binaries with e ~ 0 . 01 . 
From the integrated Euler equation, we can write this 
condition as 


Vln/i = V (In^ ) =0, (37) 

V aij 

or, by using the definitions of 70 , and 7 , 

V In = -2V In 7 . (38) 

If we decompose /3® in its inertial component /3g and its 
comoving component according to 

(3 = (3o + ^o X r + dor, (39) 

this can be written as a quadratic equation for the orbital 
angular velocity fig (neglecting the dependence of 7 on 
the orbital angular velocity fig). In practice, we solve for 


fig by projecting Eq. ( [^ 8 [ ) along the line connecting the 
center of the two starsH"^ 

The exact iterative procedure followed to solve in a 
consistent manner the constraint equations, the elliptic 
equations for (j), and the algebraic equations for h (in¬ 
cluding on-the-fly computation of fig and of the constant 
in the first integral of Euler equation) is detailed in Sec¬ 
tioning 

Once a quasi-equilibrium solution has been obtained 
by this method, lower eccentricity systems can be gen¬ 
erated by modifying fig and dg, following the methods 
developed by Pfeiffer et al. [51] . 


B. Formalism for Spinning Binaries 

We will now discuss how to alter the formalism dis¬ 
cussed above to incorporate spinning BNS. Although sev¬ 
eral formalisms have been introduced in the past Esin], 
we will follow the work of Tichy (2011)[T^. A first obvi¬ 
ous difference is that we can no longer write the velocity 
solely in terms of the gradient of a potential. Following 
Tichy, we break the velocity up into an irrotational part, 
and a new rotational part W: 

t/* =-^(a,</. + w,), (40) 

where it is natural, although not required, for W to be 
divergenceless. 

Following the assumptions stated in Tichy [15] . the con¬ 
tinuity equation becomes 


d/3*'I'^ 

Po { - + Wj) H-5i7„ -k 

' a 

= 7*4 (^d-(j) + Wi)djPo - ^ g-pg. 


(41) 


Eq. [4T] then is the same as in the irrotational case, cf. 
Eq. under the replacement di<j) —>■ di(j) + Wi. 

Taking the limit pg —0 in Eq. (41) yields the boundary 
condition at the surface of each star: 

f 4 (d.c/, + w,) djpo = (42) 


The solution of the Euler equation is no longer as 
simple as it was previously, in Eq. |20| As shown in 
Tichv('20111[16]. the solution is now 


h = Vt"- (V,(/) + W,)(V*</> + IF*), (43) 


^ Along the other directions, the enthalpy h is corrected so that 
force balance is enforced at the center of the star, according to 
the method described in m 
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where 


&+ -4a4((v,0 + 

i = - - --, (44) 


2a2 


and 


b= + + 2a^{W^cl) + w,)w\ (45) 


Finally, the method discussed previously of modifying 
the star’s angular velocity is now no longer as simple. 
The equation is modified to 


V In (a^ _ = -2V InF, 


where 

in i \r^ j cnh-y-n. / 

v'i-(? + e)(a + a) 


(46) 


(47) 


Let us now discuss the choice of the spin term, W. This 
term is, in principle, freely chosen, and so we must choose 
it so as to best represent the physical situation at hand - 
namely a uniform rotation with constant angular velocity. 
As suggested by Tichy(2011)[T6] and Tichy(2012)[T8], a 
reasonable choice for W is 


lyi ^ (48) 

where is the position vector centered at the star’s 
centre, oj^ represents an angular velocity vector and 
gijfe _ This leads to a vector field VF® with 

vanishing divergence in the conformal metric gij = Sij. 
Alternatively, one might prefer a vector field F® with 
vanishing divergence with respect to the physical met¬ 
ric gij = Owing to the conformal transformation 

properties of the divergence operator, F® is given by 

F® = 4'-®fF®. (49) 

Here, we generally use IF® as we have found that it leads 
to initial data which is closer to being in equilibrium, as 
we will further discuss in section BV El 



FIG. 1. Visualization in the x-y plane of the domain de¬ 
composition used in our initial data solve. The colour map 
represents the density of the stars. 


places all surfaces at which the solution is not smooth 
at a subdomain-boundary, which preserves the exponen¬ 
tial convergence of spectral methods. Another spherical 
shell surrounds each star. The inner shells representing 
the stars and their vicinity are embedded into a struc¬ 
ture of five concentric cylinders with three rectangular 
blocks along the axis connecting the centers of the neu¬ 
tron stars, which overlap the inner spherical shells. The 
cylinders/blocks in turn are overlapped at large radius 
by one further spherical shell centered half-way between 
the two neutron stars. Using an inverse radial mapping, 
the outer radius of the outer sphere is placed at 10^°. 

All variables are decomposed on sets of basis functions 
depending on the subdomain. The resolution of each do¬ 
main (i.e., the number of colocation points used) is cho¬ 
sen at the start of the initial data solve, and then subse¬ 
quently modified several times using an adaptive proce¬ 
dure described below. In this paper, when discussing the 
total resolution of the domain, we use the notation 

(50) 

with Ni the number of collocation points in the zth sub- 
domain. is thus the cube root of the total number 

of collaction points in all subdomains. 


C. Solving the elliptic equations 


In the previous sections, we have reduced the Ein¬ 
stein constraints, Eqs. ([T7|-([T^, as well as the conti¬ 
nuity equation (41) to elliptic equations. We solve these 


equations with the multi-domain pseudo-spectral elliptic 
solver developed in [36], as modified in [28] for matter. 
The computational domain is subdivided into individ¬ 
ual subdomains as indicated in Eig. The region near 
the center of each star is covered by a cube, overlap¬ 
ping the cube is a spherical shell which covers the outer 
layers of the star. The outer boundary of this shell is 
deformed to conform to the surface of the star. This 


D. Construction of qnasi-equilibrium initial data 


Construction of initial data for rotating binary neu¬ 
tron stars begins with selecting the physical properties 
of the system: the equation of state of nuclear matter, 
the coordinate separation d between the neutron stars, 
the baryon masses and of the two stars, and 
their spin vectors Wrotp and u>rot, 2 - We also choose the 
orbital angular frequency Hq and the initial inspiral rate 
do. 


We generally begin by setting Oq to the value for the 
orbital frequency of a similar irrotational BNS (where 
Ho is determined by the condition of quasi-circularity, 
Eq. (36)), and do = 0. These values are then adjusted 
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following the eccentricity reduction method developed by 
Pfeiffer et al. [3l]. Finally, we use a flat conformal metric, 
lij = ^ij’: and maximal slicing, K = 0. 

Once all these quantities are fixed, we need to 
solve self-consistently Eqs. 0 -([T^ for the Einstein- 
constraints, the continuity equation Eq. (41), while si¬ 
multaneously satisfying conditions to enforce the desired 
masses of the stars. To do so, we follow an iterative pro¬ 
cedure developed originally for black hole-neutron star 
binaries m- 

Eirst, we choose initial guesses for the conformal metric 
and hydrodynamical variables, using an analytical super¬ 
position of two isolated boosted neutron stars. 

We then obtain constraint-satisfying initial conditions 
by applying the following iterative procedure, where n 
represents the iteration number: 


1. Solve the nonlinear XCTS system for the set of met¬ 
ric variables X = (/3®,'k, aT), assuming fixed val¬ 
ues of the conformal source terms (E, S', J®). The 
new value of the metric variables is obtained 

from their old value X'^ and, following the relax¬ 
ation scheme used in |5S] , the solution of the XCTS 
equations X*, using 


= 0.3X* -f 0.7X”. (51) 


5. Solve the elliptic equation for the velocity potential 
4>, and obtain the next guess for cj) using the same 
relaxation method shown in Eq. |51[ 

6 . Check whether all equations are satisfied to the de¬ 
sired accuracy. If yes, proceed. If no, return to 
Step[2 

7. Compute the truncation error of the current so¬ 
lution by examining the spectral expansion of the 
XCTS variables. If this truncation error is un¬ 
desirably large (typically, if it is > 10“®), then 
adjust the number of grid-points in the domain- 
decomposition and return to Step The adjust¬ 
ment is based on the desired target truncation eror 
and the measured convergence rate of the solution, 
cf. ng. 


E. Quasi-Local Angular Momentum 

The goal of the present paper is to construct spinning 
BNS initial data and to evolve it. Therefore, we need 
diagnostics to measure the NS spin, for which we use 
techniques originally developed for black holes. It is com¬ 
mon to discuss the spins of black holes in terms of their 
dimensionless spin x, 


2. Locate the surface of each star. Representing 
the surface in polar coordinates centered on each 
star as i?"(0, 0), we determine i?” to satisfy [55] 
h{R'2{0, (j)), 9, (j)) = 1. To ensures that the grid¬ 
boundary Rb converges to the surface of the star, 
we occassionally modify the numerical grid such 
that Rb{9,(j)) = R^{9,(j)). Because this requires 
a re-initialization of the elliptic solver, the grid is 
only modified if the stellar surface has settled down, 
specifically, if 

\\R:-R:-^\\<o.i\\R:-Rb\\. ( 52 ) 


X = 


s 


(54) 


Here, S is the angular momentum of the black hole, and 
M is its Christodoulou mass [5^ . 


= Ml 


4M2/ 


(55) 


The irreducible mass is defined base d on th e area 
of the hole’s apparent horizon, Mi„ = y/AJlGTr. The 
angular momentum is computed with a surface integral 
over the apparent horizon 


Here || . II 2 denotes the L2-norm over the surface. 

3. Eor each neutron star, fix the constant in Euler’s 
first integral so that the baryon mass of the neutron 
star matches the desired value. We compute the 
baryon mass as a function of the Euler constant C 
through 

and utilize the secant method to drive the mass to 
the desired value. 


4. If desired, adjust the orbital frequency to ensure 
force-balance at the center of each star by solv¬ 
ing Eq. (38). This step is skipped if the orbital 
frequency is fixed through iterative eccentricity re¬ 
moval, cf. Sec. IV B[ 


S = — (t (j)^s^K,jdA 

OTT 7^ 


(56) 


where R is the black hole’s apparent horizon, is the 
outward-pointing unit-normal to R within the t = const 
hypersurface, and is an azimuthal vector field tangent 
to R. Eor spacetimes with axisymmetry, c/i® should be 
chosen as the rotational Killing vector. In spacetimes 
without an exact rotational symmetry (e.g. the space- 
time of a binary black hole system), one substitutes an 
approximate Killing vector [321 011 (AKV). Ref. [33] in¬ 
troduces a minimization principle to define <()*, resulting 
in an Eigenvalue problem. The three eigenvectors with 
the lowest eigenvalues (i.e. smallest shear) are taken and 
used to compute the three components of the spin. 

In this paper, we explore the application of quasi-local 
spin measures to neutron stars. In the absence of appar¬ 
ent horizons 77, we need to choose different surface(s) to 
evaluate Eq. (56). 
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When constructing initial data, the stellar surface S is 
already determined, so one obvious choice is to integrate 
over the stellar surface S. To estimate the ambiguity 
in quasi-local spin, we furthermore compute S by inte¬ 
grating over coordinate spheres with radii ranging from 
just outside S to larger by about 70%. During the evo¬ 
lution, the stars change shape and may even loose mass 
in tidal tails. Because of these complications, the SpEC 
evolution code does not track the location of the stellar 
surface during the evolution, and we shall only monitor 
S on coordinate spheres. 

It is useful to compute a dimensionless spin y, for in¬ 
stance, for post-Newtonian comparisons. In the absence 
of a horizon, Eq. (55) is meaningless and we need a differ¬ 


ent choice for the mass-normalization. Instead, we nor¬ 
malize by each star’s ADM mass, Madm, i.e. 


X = 


s 

/\/f2 


(57) 


The ADM mass is determined by computing the ADM 
mass of an equilibrium configuration of a single uniformly 
rotating polytrope in isolation with the same baryon 
mass and angular momentum as those measured in our 
binary systems. 

The results of the quasi-local spin measures are de¬ 
scribed in section IIID[ which shows that this procedure 
is numerically robust. 

Finally, let us discuss, from an order of magnitude per¬ 
spective, how the star’s dimensionless spin is related to its 
more commonly used physical properties. We start with 
the Newtonian relation S = 2ttI/P between angular mo¬ 
mentum S, moment of inertia /, and rotational period 
P. Writing further I = f M, with the dimensionless 
constant / depending on the stellar density profile, we 
have 


27rc fR^ 

^ ~g¥m 

= 0.48 


(^ ] 

( ^ ] 

2 ^ 

r ^ ] 


r ^ ^ 

V0.33/ 

V 12km J 

KiaMqJ 

Urns/ 


(58) 


The factor c/G arises from the transition to geometric 
units. 

This -quite simplistic- estimate shows that millisecond 
pulsars will have appreciable dimensionless spin y. Cen¬ 
trifugal breakup of rapidly rotating neutron stars hap¬ 
pens at a dimensionless spin in the range 0.65 — 0.70 [iH] . 
with only small dependence on the equation of state and 
neutron star mass. Ansorg et al [15] studied in detail 
r = 2 polytropes, the equation of state we use here. 
Ref. [46| finds a dimensionless spin at mass-shedding of 
X = 0.57. 


d 

I 

C 

o 

U 


:3 

W 



_,_^_,_I_,_I_^_L_ 

0 20 40 60 80 

iteration # 


FIG. 2. Convergence of the Euler constant during iteration 
at the lowest resolution RO. The inset shows the difference 
between values at subsequent iterations. 


for BNS systems with arbitrary spins. As discussed in 
section m our code consists of a solver that runs for a 
number of iterations at constant resolution, and then the 
resolution is increased and this process restarts. We will 
therefore demonstrate that appropriate quantities con¬ 
verge with both the iterations of iterative scheme de¬ 
scribed above in Section fll PI and with resolution as the 
resolution increases. 


A. Convergence of the Iterative Procedure 

At each step of the iterative procedure, the Euler con¬ 
stant of each star is modified to achieve a desired stellar 
baryon mass, based on the current matter distribution 
inside the star. We expect that the Euler constant con¬ 
verges during the iterations at a fixed resolution. Fig¬ 
ure [^ shows the behavior of the Euler Constant during 
iterations at the lowest initial data resolution, RO. We 
show three runs of interest, one with large aligned spins 
(S.4z), one with large precessing spin (S.4x), and one 
with small anti-aligned spins (S-.05z). The properties 
of these configurations are shown in table |T| In all three 
cases we see agreement between neighboring iterations 
at the 10 “® — 10 “® level by the end of iterating at this 
resolution. At the highest resolutions, these differences 
are down to, typically, the 10“® — 10“^° level. This can 
be compared to Fig. 3 of m- Although not shown here, 
other free quantities converge similarly to the Euler con¬ 
stant. 


B. Convergence of the Solution 

III. INITIAL DATA RESULTS 

Having established that our iterative procedure con- 
In this section, we will demonstrate that our code verges as intended, we now turn out attention to the 

can robustly construct contraint-satisfying initial data convergence of the solution with resolution. To demon- 


























Name 

Mns 

u; 

Do 

Do X 10® 

ho X 10® 

Madm 

X 

S.4z 

1.7745 

0.01525Z 

47.2 

5.09594 

-1.75 

1.648 

0.37652 

S-.05Z 

1.7745 

-0.00273£ 

47.2 

5.11769 

-1.71 

1.640 

-0.050182 

S.4x 

1.7745 

0.01525* 

47.2 

5.10064 

-2.36 

1.648 

0.3714* 


TABLE I. Parameters for the initial data sets used in testing the initial data solver: M^s and cu® are baryon mass and 
rotational parameter for either neutron star (the same values are used); Do, flo and ho represent coordinate separation between 
the centers of the stars, the orbital frequency, and the radial expansion; x is the dimensionless spin vector computed from the 
initial data set. In each case we use a polytropic equation of state, P = Kpo, with P = 2 and k = 123.6. 


strate it, we will look at the Hamiltonian and momentum 
constraints, and the differences between measured phys¬ 
ical quantities - the ADM energy and ADM angular mo¬ 
mentum, and the surface fitting coefficients of the stars. 
As our initial data representation is fully spectral, we ex¬ 
pect that these quantities should converge exponentially 
with resolution. Note that when we discuss the value of 
a quantity at a certain resolution, we are referring to the 
value of that quantity after the final iterative step at that 
resolution. 

Figure shows the convergence of the Hamiltonian 
constraint and the Momentum constraint for our three 
runs of interest. These are computed during the last it¬ 
erative solve at each resolution. The data plotted are 
computed as 


1 II 

(59) 

Rp II 

(60) 


Here R^j and Rp denote the residuals of Eqs. (18) 
and respectively, and 11.11 represents the root-mean- 
square value over grid-points of the entire computational 
grid. This plot demonstrates that our initial data solver 
converges exponentially with resolution, even for very 
high spins, which gives confidence that we are indeed 
correctly solving the Einstein Field Equations. 

The surface of the star is represented by a spherical 
harmonic expansion: 


‘■max •)'' ‘'max 


Rs{0,(j))= ^ ClmYi^{9,(j}), 


(61) 


l,m 


where Zmax = Wmax = 11, unless stated otherwise. The 
stellar surface is located by finding a constant enthalpy 
surface, cf. Sec. H D[ and the spherical subdomains that 
cover the star are deformed to conform to Rs{0,(j)). To 
establish convergence of the position of the stellar surface 
we introduce the quantity 


Ac{i) = 


l{l + l)\ 


‘'max 1 '' ‘■max 


(62) 


Here i refers to the resolution in the initial data, and 
N refers to the final resolution. Figure]^ plots Ac(i) vs. 



FIG. 3. Hamiltonian and Momentum constraints as a function 
of resolution N. We see exponential convergence in all cases. 


O 



FIG. 4. Gonvergence of the lo cat ion of the stellar surface. 
Plotted is Ac as dehned in Eq.( 62 I, for three representative 
conhgurations. 


resolution. The surface location converges exponentially 
to better than 10“®. 

Finally, we assess the overall convergence of the so¬ 
lution through the global quantities Aadm and 
The surface integrals at infinity in these two quantities 
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FIG. 5. Convergence of ADM-energy and the magnitude 
of the ADM-angular momentum. Shown are the fractional 
differences between neighboring resolutions, as a function of 
the lower resolution. 


are recast using Gauss’ law (cf. [^l: 


£^adm ~ ^ dS. 


2tt 


= -— i 5)d,-^ds^ + — / (63) 


27r 


and 


•^ADM = ^ / {^Ky^ - yK^^) dSj 


= — ^ {xKy, - yK^i) 5^^ dSj. (64) 


Here V is the volume outside S, and the integrals are eval¬ 
uated in the flat conformal space. To obtain the other 
components of Jadmj cyclically permute the indices x,y,z. 
We define the quantities AE and A J as the absolute frac¬ 
tional difference in these quantities between the current 
resolution and the next highest resolution. These are 
plotted in figure In general, we find agreement at the 
10 “^ — 10“® level by the final resolution. 


C. Convergence of the quasi-local spin 

We now turn to the angular momentum of the neutron 
stars, as measured with quasi-local angular momentum 
integrals on the stellar surface. We will discuss dimen¬ 
sionless spins X, which depend on two distinct numeri¬ 
cal resolutions: First, the resolution of the 3-dimensional 
grid used for solving the initial value equations. This 



FIG. 6. Convergence of the quasi-local spin computation. 
Top panel: difference of spin computed at resolution N with 
the spin computed at the highest resolution. Bottom panel: 
Difference between spins computed at different resolution L 
of the spin-computation. For S-.5z, we achieve an accuracy 
of ~ 10“^, whereas for S.4z and S.4x, the accuracy is ~ 10“"^ 
due to finite L. 


resolution is specified in terms of N, the total number 
of grid-points. Second, the resolution used when solv¬ 
ing the eigenvalue problem for approximate Killing vec¬ 
tors on the 2-dimensional surface, as given by L, the 
expansion order in spherical harmonics of the surface- 
parameterization rs{B, (j)) = J2f=o (j)). 

Throughout this paper, we use L = 11. The top panel 
of Fig. shows convergence of x with grid-resolution N, 
at fixed L = 11. We find near exponential convergence. 

The influence of our choice L = 11 is examined in the 
lower panel by computing the quasi-local spin at lower 
resolution L = 8 and at higher resolution L = 14. Chang¬ 
ing L impacts x by 10“® for the low-spin case S-. 05z, 
and by ~ 10“"^ for the high-spin cases S. 4z and S. 4x. For 
the high-spin cases, the spin measurement is convergent 
with increasing L, and the finite value of L dominates the 
error budget. For the low-spin case, numerical truncation 
error dominates the error budget and convergence with 
L is not visible. High NS spin leads to a more distorted 
stellar surface, and so a fixed L = 11 yields a spin result 
of lower accuracy. However, in all cases the the numeri¬ 
cal errors of our spin measurements are still neglible for 
our purposes. 


D. Quasi-local Spin 

As discussed in section El we use a quasi-local spin 
to define the angular momentum carried by each neutron 
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FIG. 7. Stellar cross-sections in the X-Z plane for a series 
of different spins, aligned with the z axis, demonstrating that 
they bulge at the equator in the expected way with increasing 
spin. 



FIG. 8. Dimensionless angular momentum x S'® ^ function of 
fl for a series of spin-aligned initial data sets with the same 
physical parameters as our runs of interest. We see, as ex¬ 
pected, a linear relation between x and at low-spins, which 
eventually becomes non-linear at higher spins. 


star. To our knowledge, this is the first application of this 
method to neutron stars in binaries. 

In this section, we explore properties of the rotating 
BNS initial datasets and the employed quasi-local spin 
diagnostic. 

To explore the spin-dependence of BNS initial data 
sets, we construct a sequence of equal-mass, equal-spin 
BNS binaries, with spins parallel to the orbital angular 
momentum. We fix the initial data parameters M^g, Dq, 
Ho and do to their values for a configuration that we will 
also evolve below (specifically, S.4z - Eccl) 

Figure shows cross-sections through one of the neu¬ 


tron stars in the xz-plane, i.e. a plane orthogonal to the 
orbital plane which is intersecting the centers of both 
stars. With increasing spin, the stars develop an increas¬ 
ing equatorial bulge, an expected consequence of centrifu¬ 
gal forces. 

Figure [^presents the dimensionless spin of either neu¬ 
tron star as a function of w. x increases monotonically 
with the rotation parameter w. The spin x increases lin¬ 
early with w for small w. For larger w, the dependence 
steepens, as the increasing equatorial radius of the stars 
increase the moment of inertia [15]. 

For uj = 0.01625Mq^ we achieve x = 0.432, the largest 
spin we are able to construct. This is reasonably close 
to the theoretical maximum value for T = 2 polytropes, 
X 0.57 [15]. Above uj = 0.01625Mq^, the initial data 
code fails to converge. The steepening of the x vs. uj 
curve is reminiscent of features related to non-uniqueness 
of solutions of the extended conformal thin sandwich 
equations miiMT], and it is possible that our failure 
to find solutions originates in an analogous break-down 
of the uniqueness of solutions of the constraint equations. 

While the focus of our investigation lies on rotating NS, 
we note that for a; = 0 our data-sets reduce to the stan¬ 
dard formalism for irrotational NS. For w = 0, we find a 
quasi-local spin of the neutron stars is x = 2 x 10“^. This 
is the first rigorous measurement of the residual spin of ir¬ 
rotational BNS. Residual spin is, for instance, important 
for the construction and validation of waveform models 
for compact object binaries. The analysis in Ref. [^ in¬ 
dicates that spins of order 10“^ lead to a dephasing of 
about O.Olradians during the last dozen of inspiral orbits. 
This value is significantly smaller than the phase accu¬ 
racy obtained by current BNS simulations, and so the 
residual spin is presently not a limiting factor for studies 
like [53l - f55] . 

Finally, we demonstrate that the surface on which we 
compute the quasi-local spin, does not significantly im¬ 
pact the spin we measure: We choose coordinate spheres 
centered on the neutron star with radius R, and compute 
the quasi-local spin using these surfaces, rather than the 
stellar surface. 

In Fig. we plot the spin measured on various R = 
const surfaces, for three different values of w, from the 
same sequences shown in Fig.[^ 

The circles denote spins extracted on coordinate 
spheres. The asterisks indicate the spins computed on 
the stellar surface. The asterisk is plotted at R — i?oq, 
the equatorial radius of the neutron star under considera¬ 
tion. We find good agreement between spins extracted on 
coordinate spheres and the spin extracted on the stellar 
surface, as long as R > i?eq- The maximum disagreement 
is seen in the high spin curve, where the two spins differ 
by ~ 10-2. 

For R < Rcq, the coordinate extraction sphere inter¬ 
sects the outer layers of the neutron star and no longer 
encompasses the entire matter and angular momentum 
of the star. Therefore, x(^) shows a pronounced de¬ 
cline for R < Rcq for each of the three initial-data sets 









11 



FIG. 9. Dimensionless spin x measured on coordinate spheres 
with radius R for three different aligned spin BNS systems. 
The asterik denotes the spin measured on the (non-spherical) 
stellar surface. Circles to the right of the asterik represent co¬ 
ordinate spheres entirely outside the neutron star, and circles 
on the left of the asterik indicate spin measurement surfaces 
that intersect the star or are entirely located inside the star. 


considered in Fig. ^ For R > i?oq, x(i?) continues 
to increase slightly, for instance, for the middle curve, 
x(i? = 9) = 0.202 whereas x(i? = 11) = 0.204. 

In summary. Fig. shows that the quasi-local spin ex¬ 
tracted on coordinate spheres can serve as a good ap¬ 
proximation of the quasi-local spin extracted on the stel¬ 
lar surface (as long as the coordiate sphere is outside the 
star, of course). 

This is important because during evolutions of the bi¬ 
nary, we do not track the surface of the star. Instead, we 
will compute the spin on coordinate spheres, similarly to 
Fig.0 


IV. EVOLUTION RESULTS 

We now evolve the three configurations discussed in 
Sec. |III| As indicated in Table all three configurations 
are equal-mass binaries, with individual ADM masses M* 
(in isolation) of 1.64M0 or 1.648M0 at initial separation 
oi D = 47.2M0, and using a polytropic equation of state 
with F = 2.0 and k = 123.6. Both stars have equal 
spins, and the three configurations differ in spin mag¬ 
nitude and spin direction. Configuration S-0.05z has 
spin-magnitudes ~ 0.05 anti-aligned with the orbital an¬ 
gular momentum, and the confiurations S. 4z and S. 4x 
have spin magnitudes near 0.4, along the z-axis and x- 
axis, respectively. 

Each configuration is evolved through > 10 orbits, into 
the late-inspiral. In this paper we focus on the inspiral 
of the neutron stars. Table |TT] summarizes parameters for 
these runs. 


Name 

k 

e 

X 

MHz) 

Norb 

tf (ms) 

S.4z 

0 ,1,2 

< O.OOl 

0.3815 

167.7 

11.8 

56.0 

S-.05Z 

0 ,1,2 

0.0006 

-0.0505 

165.4 

12.5 

56.3 

S.4x 

0,1 

< 0.002 

0.375X 

164.8 

9.1 

45.7 


TABLE II. Information about our three evolutions, k indi¬ 
cates the numerical resolutions on which a simulation is per¬ 
formed, e indicates the smallest achieved orbital eccentricity. 
X and /o are the dimensionless spins at t = 0 and the initial 
orbital frequency. Finally, Aorb and tf represent the number 
of orbits the configuration was evolved for, and the evolution 
time. 

A. Evolution Code 

In our evolution code, SpEC we use a mixed 

spectral - finite-difference approach to solving the Ein¬ 
stein Field Equations coupled to general relativistic hy¬ 
drodynamics equations. The equations for the space- 
time metric, are solved on a spectral grid, while the 
fluid equations are solved on a finite difference grid, us¬ 
ing a high-resolution shock-capturing scheme. We use 
a WENO [5S1 in?] reconstruction method to reconstruct 
primitive variables, and an HLL Riemann solver [68] to 
compute numerical fluxes at interfaces. Integration is 
done using a 3rd order Runge-Kutta method with an 
adaptive stepsize. We interpolate between the hydro and 
spectral grids at the end of each full time step, interpo¬ 
lating in time to provide data during the Runge-Kutta 
substeps (see [smsMi] for a more detailed description 
of the method). 

Each star is contained in a separate cubical finite dif¬ 
ference grid that does not overlap with that of the other 
star. The sides of the grids are initially 1.25 times 
the stars’ diameters. We use grids that contain 97^, 
123^ and 155^ points for resolutions k = 0,1,2, re- 
spectiveljj^ These resolutions correspond to linear grid¬ 
spacing of 340 m, 268 m and 213 m respectively for the 
S. 4z case. The precessing evolution S. 4x uses similar 
grid-spacing, whereas the anti-aligned run S-. 05z has a 
slightly smaller grid-spacing because the stars themselves 
are smaller. The region outside the NS but inside the fi¬ 
nite difference grid is filled with a low density atmosphere 
with p — 10“^^Mq^. The motion of the NSs is monitored 
by computing the centroids of the NS mass distributions 

^CM = J x"u°Pay/- (65) 

for each of the grid patches containing a NS. 

The grids are rotated and their separation rescaled 
to keep the centers of the NS at constant grid- 
coordinates [SIl [Ml Eg. As the physical separa- 


^ For aligned-spin configurations S-.05z and S.4z, we take advan¬ 
tage of, and enforce, z-symmetry, which halves the number of 
grid-points along the z-axis. 
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tion between the stars decreases, the rescaling of grid- 
coordinates therefore causes the size of the stars to in¬ 
crease in grid-coordinates. In order to avoid the stellar 
surfaces expanding beyond the geometric size of the finite 
difference grid, we monitor the matter flux leaving this 
grid along the x, y, and z-direction. If the matter flux 
is too large along a certain axis, we expand the grid in 
that direction. This procedure allows us to dynamically 
choose the optimal grid-size that limits matter loss to a 
small, user-specified level. When changing the size of the 
hydro grid, the number of grid-points is kept constant, 
so this process changes the effective resolution during the 
evolution. 

The Einstein field equations are solved on a spectral 
grid using basis-functions appropriate for the shape of 
each subdomain. For rectangular blocks, Cheybyshev 
polynomials are used along each axis; for a spherical 
shell (i.e. where the center is excised), spherical har¬ 
monics in angles, and Chebyshev polynomial in radius 
are employed; and for an open cylinder (i.e. with the re¬ 
gion near the axis excised), Chebyshev polynomials and 
a Fourier series. For full spheres and filled cylinders, 
multi-dimensional basis-functions respecting the continu¬ 
ity conditions at the orign/axis are employed [731 [Ti]. For 
more details see m- 

More specifically, our spectral grid, the central region 
of each star is covered by a filled sphere located at the 
center of the star. These have spherical harmonic modes 
up to L = 12 -I- 2k. The radial basis-functions are one¬ 
sided Jacobi polynomials with 7 -|- fc collocation points. 
The filled spheres are surrounded by eight other spherical 
shells with the same radial and angular resolutions. At 
the start of the evolution, the stellar surface is generally 
located inside the third shell. The far field region is cov¬ 
ered by 20 spherical shells starting at 1.5 times the inital 
binary separation and going out to 40 times that separa¬ 
tion. These shells have angular resolution L — 9 + 2k and 
radial resolution 6-|-fc. The region between the innermost 
shell and the stars is covered by a set of cylindrical shells 
and filled cylinders. 

We use a generalized harmonic evolution system inni 
EillTS] with coordinates such that they satisfy a wave 
equation 

( 66 ) 

for some freely-specifiable source function The ini¬ 
tial source function is determined by the intial 

data, assuming that the time derivatives of the lapse and 
shift functions initially vanish in the corotating frame. 
We then transition to a pure harmonic gauge, = 0 
by using a transition function, i.e. 

( 67 ) 

The timescale r is determined by r = 2 a/#/( 2M*). This 
is slow enough to avoid numerical gauge artifacts in the 
simulations. 



FIG. 10. The binary separation as a function 

of time. Shown are three eccentricity removal iterations 
(Eccl,Ecc2,Ecc3) for each of the three configurations stud¬ 
ied. The data for S-. 05z and S. 4z is offset vertically by 6 
and 3, respectively, for clarity of plotting. 


B. Eccentricity Removal 

Gravitational wave emission reduces orbital eccentric¬ 
ity rapidly during a GW-driven inspiral [131E] . There¬ 
fore inspiraling binary neutron stars are expected to have 
essentially vanishing orbital eccentricity in their late in¬ 
spiral, unless they recently underwent dynamical inter¬ 
actions. Our goal is to model non-eccentric inspirals. In 
this subsection we demonstrate that we can indeed con¬ 
trol and reduce orbital eccentricty, using the techniques 
developed for BH-BH binaries [31E1I77] and also ap¬ 
plied to BH-NS binaries [35] . 

For fixed binary parameters (masses, spins), and fixed 
initial separation Dq, the initial orbit of the binary is 
determined by two remaining parameters: The initial or¬ 
bital frequency Og, and the initial radial velocity, which 
we describe through an expansion parameter dg = r/r. 
These two parameters will encode orbital eccentricty and 
phase of periastron, and our goal is to determine these 
parameters to reduce orbital eccentricity. We accomplish 
this using an iterative procedure first introduced for bi¬ 
nary black holes |52l [TT] . An initial data set is evolved 
for a few orbits, the resulting orbital dynamics are an¬ 
alyzed, and then the initial data parameters Ug and dg 
are adjusted. 

For binary neutron stars, we initialize the first itera¬ 
tion of eccentricity removal, with dg = 0 and use Ug de¬ 
termined from irrotational BNS initial data, based on the 
equilibrium condition in Eq. |38| Evolutions with these 
choices are labeled with the suffix “Fed”, and show not- 
icable variations in the separation between the two NS, 
cf. the solid black lines in Fig. 

We compute the trajectories of the centers of mass of 
each star, as determined by Eq. Ci(t) and C 2 (t), and 
using the relative separation f = C 2 (t) — cj(t), compute 
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Name 

n X 10® 

do X 10® 

e 

S.4z - Eccl 

5.10538 

0 

0.006 

S.4z - Ecc2 

5.09591 

-1.60 

< 0.001 

S.4z - Ecc3 

5.09594 

-1.75 

< 0.001 

S-.05Z- Eccl 

5.10538 

0 

0.008 

S-.05Z- Ecc2 

5.11561 

0 

0.004 

S-.05Z - Ecc3 

5.11769 

-1.71 

0.0006 

S.4x - Eccl 

5.10538 

0 

0.007 

S . 4x - Ecc2 

5.10429 

-2.27 

0.004 

S . 4x - Ecc3 

5.10064 

-2.36 

< 0.002 


TABLE III. Eccentricity removal for the three main runs 
discussed in this paper. Only initial orbital frequency flo 
and initial radial expansion factor ao are changed between 
different EccN iterations. Recall that these quantities have 
units of 



Time (IVy 



Time (M^ 


FIG. 11. The derivative of the binary orbital frequency as a 


function of time for different levels of eccentriccity reduction 


for our three runs of interest. Note that dQ/dt has units of 


the orbital frequency 


n{t) 


W) X f(t)\ 

r{tY 


( 68 ) 


where an over-dot indicates a numerical time-derivative. 
Finally, we compute f2(t) and fit it to a function of the 
form 


(lit) =Ai{tc-t) -I-^ 2(4 - t) 

-l- Bq cos {Bit + ^ 3 )- (69) 


The power law parts of this fit represent the orbital de¬ 
cay due to the emission of gravitational waves, while the 
oscillatory part represents the eccentric part of the orbit. 
We then update flo and oq with the formulae (see [77] 


FIG. 12. Convergence of li(t). Shown are li(t) at three dif¬ 
ferent numerical resolutions (fc = 0,1, 2) for the final, lowest- 
eccentricity initial data. The oscillations in fl{t) are evidently 
not caused by numerical truncation error. Note that Omega 
has units of 


for a detailed overview) 


0 , 0 BqBi , ^ 

iZo ilo ^^2 

(70) 

• ^ • _L R 

+ 7 ^7^ COS R 3 . 

Z\ Zq 

(71) 


We repeat this procedure twice, resulting in simulations 
with suffix Ecc2 and Ecc3. Table imi summarizes the 
orbital parameters for the individual simulations, and 
Figs. [To] and 11 illustrate the efficacy of the procedure 
through plots of separation and time-derivative of orbital 
frequency. The eccentricity is successfully reducued from 
e ^ 1% to ^ 0.1%. After two eccentricity reduction it¬ 
erations, variations in n{t) are so small that they are 
no longer disce rnib le from higher-frequency oscillations 
in n{t), cf. Fig. 
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The high freuqency oscillations in il{t) are caused 
by the quasi-normal ringing of the neutron stars, as 
discussed in detail below in Sec. |IVE| Here, we only 
note that these oscillations are convergently resolved, cf. 
Fig. |12[ and are therefore a genuine feature of our ini¬ 
tial data. Figure also confirms that the lowest reso¬ 
lution (fc = 0) gives adequate resolution for eccentricity 
removal. 


The eccentricity removal algorithm attempts to isolate 
variations on the orbital time-scale as the signature of 
eccentricity. For S.4z - Ecc2, it reports e = 0.0005 and 
for S.4z - Ecc3, e = 0.0002. However, given the large 
amplitude of the QN mode ringing, we consider these 
estimates unreliable, and therefore quote an upper bound 
of 0.001 in Table jlllj Similarly, for S. 4x - Ecc3, the fitting 
reports e = 0.001, and we quote a conservative upper 
bound of 0.002. 
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FIG. 13. The spin measured on multiple coordinate spheres 
for the S. 4z run. 


C. Aligned spin BNS evolutions: NS Spin 


In this section, we will discuss the measurement of 
spins during our evolutions for the non-precessing cases, 
S.4z and S-.05z. Aligned spin binaries do not pre- 
cess. Combined with the low viscosity we expect the 
NS spins to stay approximately constant during the evo¬ 
lutions. These systems therefore serve as a test on our 
spin diagnostics during the evolutions. In this section, 
and through the rest of this paper, we always use the 
final eccentricity reduction, “Ecc3”. For brevity, we will 
omit the suffix “-Ecc3”, and refer to the runs simply as 
S-.OSz, etc. 


We do not track the surface of the star during the evo¬ 
lution. Instead we simply evaluate the quasi-local spin 
of the stars on coordinate spheres in the frame comov¬ 
ing with the binary. We must therefore verify that the 
spin measured is largely independent of the radius of the 
sphere, and that it is maintained during the evolutions 
at the value consistent with that in the initial data. Fig¬ 
ure established that coordinate spheres can be used to 
extract the quasi-local spin in the initial data. Figure 
shows the results for the high-spin simulation S. 4x dur¬ 
ing the inspiral. 

For coordinate spheres with radii R = 11.28M0 to i? = 
16.93M0 in grid coordinates, the spins remain roughly 
constant in time. The different extraction spheres yield 
spins that agree to about 1%, with a consistent trend that 
larger extraction spheres result in slightly larger spins 
(as already observed in the initial data). The horizontal 
dashed line in Fig. [T^ indicates the spin measured on the 
stellar surface (i.e. not on a coordinate sphere) in the 
initial data. We thus find very good agreement between 
all spin measurements, and conclude that the quasi-local 
spin is reliable to about 1%. 


The extraction sphere R = 9.87Mq in Fig. 13 



Time (M^ 


FIG. 14. Neutron star spin during the two aligned-spin 
evolutions. Shown are three different numerical resolutions, 
k = 0 (lowest), k = 1, and fc = 2 (highest). The asterisk 
indicates the spin measured on the stellar surface in the initial 
data. 


tersects the outer layers of the neutron star. Because 
the quasi-local spin captures only the angular momen¬ 
tum within the extraction sphere, the value measured on 
R = 9.87Mq falls as our comoving grid-coordinates cause 
this coordinate sphere to slowly move deeper into the in¬ 
terior of the star. This behavior, again, is consistent with 
Fig.i 

These tests of using multiple coordinate spheres were 
only run for about half of the inspiral - enough to es¬ 
tablish that the method is robust. Subsequently, we re¬ 
port spins measured on the largest coordinate sphere, 
i?= 16.93M©. 


The full behavior of the spin during the inspiral is 
shown in figurefor both the S.4z and S-.05z runs. 
Comparing the spin at different resolutions, we note that 
the data for fc = 1 and fc = 2 are much closer to each other 
than compared to fc = 0, indicating numerical conver¬ 
gence. We note that the impact of numerical resolution 
(as shown in Fig. 14) is small compared to the uncertainty 
inherent from the choice of extraction sphere, cf. Fig.[T^ 
We also note that for the first lOOOOM© of the run, the 
measured spin behaves as a constant, as expected, albeit 
with some small oscillations. However, afterward, we no¬ 
tice the absolute value of the spin starts to decrease in 
both cases. Finally, we note that in both cases, the spin 
measured in the inital data on the stellar surface is within 
Ay = 0.008 of the spin measured during the evolution. 

Finally, we compute the orbital phase 




n(t') dt', 


(72) 
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FIG. 15. Accumulated orbital phase as a function of time for 
our anti-aligned, S-.05z, and aligned, S.4z, runs. The dashed 
lines are Taylor T4 PN simulations. The PN simulations were 
matched to NR in the intervals [1109,3956] and [2090,4904] 
respectively. Qualitatively, there is excellent agreement with 
the numerical data. The lower panel shows the difference 
A0(t) = 0NR(t) — (^PN(t). 


where the orbital frequency Q(t) is given by Eq. ( 681 . 
The result is plotted in Fig. along with the 

Post-Newtonian prediction for the same binary param¬ 
eters (spins, masses and initial orbital frequencies). We 
use the Taylor T4 model (see e.g.,[^) at 3.5PN order ex¬ 
pansion, with no tidal terms added, using the matching 
techniques described in m- We find excellent qualita¬ 
tive agreement in both cases, thereby giving additional 
evidence that our numerical simulations are working as 
expected. We do hnd large late time growth in the phase 
difference, however this is expected because we do not 
model tidal effects, which become increasingly important 
at late times, in our Post-Newtonian equations. 

Figure shows the gravitational waveforms for our 
two non-precessing simulations. We extract the waves 
on a sphere of radius R — 627Mq. 



FIG. 16. The gravitational waveforms for our anti-aligned, 
S-.05z, and aligned, S.4z runs. The black curve represents 
the real part of the waveform, 5R(/i2,2) while the orange curve 
represents the magnitude of the waveform. 



FIG. 17. Spin-components of one of the neutron stars during 
the precessing simulation (thick, solid lines). The dotted and 
dashed lines represent the unmatched and matched PN results 
respectively. The agreement between PN and NR is good for 
both PN simulations. The orbital frequency was evolved using 
the Taylor T4 approximant. The matching was done in the 
interval [1892,4575]. 


D. Precession 

We now turn to the precessing simulation, S.4x. Fig¬ 
ure shows the components of the spin-vector x of one 
of the neutron stars, as a function of time. The quasi¬ 
local spin diagnostic returns a spin with nearly constant 
magnitude, varying only by ±0.002 around its average 
value 0.370. The spin components clearly precess, with 
the dominant motion in the xy-plane (the initial orbital 
plane), with the simulation completing about 2/3 of a 
precession cycle. A z-component of the NS spin also ap¬ 


pears, indicating precession of the neutron star spin out 
of the initial orbital plane. 

Fig. |17| shows a comparison of spin precession between 
numerical relativity and Post-Newtonain theory. We per¬ 
form this comparison using the matching tecnhique in 
m- This gives very good agreement between PN (dot¬ 
ted) and NR (solid) as shown by Fig. 17 The NS spins in¬ 
deed precess as expected, thus confirming both the qual¬ 
ity of quasi-local spin measures, as well as the perfor¬ 
mance of the PN equations. Note that z-component of 
the spin in the NR data undergoes oscillations that are 
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FIG. 18. Components of the orbital frequency vector fl. 
Thick solid lines represent the processing BNS simulation and 
thin dashed lines represent the matched post-Newtonian sim¬ 
ulations. The inclination reaches 5 = 0.34rad at t = 7600Mq. 




Time (IVy 


FIG. 19. Gravitational waveforms of our precessing run. 
Shown are the {I, m) = (2,2) and (2,1) modes, as extracted 
in a spherical harmonic decomposition aligned with the z-axis. 
The emergence of the (2,1) mode indicates precession of the 
orbital plane away from the xy-plane. 


unmodelled by PN. These occur on a timescale of half the 
orbital timescale. Similar effects were found in [^. The 
origin of these oscillations remains unclear. The preces¬ 


sion of the orbital angular frequency is shown in Fig. 18 


We find substantial precession away from the initial di¬ 
rection of the orbital frequency flo oc z, with the angle 
5 between 17 (t) and the z-axis reaching 20°. Once again, 
the PN equations reproduce the precession features suc¬ 
cessfully. 

Finally, Fig. 19 shows the (2,2) and the (2,1) spheri¬ 
cal harmonic modes of the gravitational wave-strain ex¬ 
tracted at an extraction surface of radius R = GdTM©. 
The (l,m)=(2,l) mode would be identically zero for an 
equal-mass aligned spin binary with orbital frequency 


FIG. 20. The maximum density p{t) in each of our runs, 
normalized by the initial maximum density p(0). The inset 
shows an enlargement of all three runs, illustrating that the 
oscillations are more pronounced in the high-spin simulations. 


parallel to the z-axis, so the emergence of this mode once 
again indicates precession in this binary. 


E. Stellar Oscillations 


The rotating neutron stars constructed here show oscil¬ 
lations in the central density, as plotted in Fig.[^ In the 
low spin run, the density oscillations have a peak-to-peak 
amplitude of about 0.6%, whereas in the high-spin runs 
(S.4z and S.4x), the density oscillations reach a peak- 
to-peak amplitude of 20%. The two high-spin simula¬ 
tions show oscillations of nearly the same amplitude and 
frequency, therefore oscillating nearly in phase through¬ 
out the entire inspiral. The oscillation-period is about 
177Mq ~ 0.87ms, i.e. giving a frequency of 1.15kHz. 
It remains constant throughout the inspiral. The low- 
spin run S-0.5z exhibits a slightly smaller oscillation pe¬ 
riod of about P « 170Mq « 0.84ms, i.e. a frequency of 
« 1.19kHz. 

To investigate the spectrum of the density oscillations, 
we perform a Fourier-transform on p{t). The result is 
shown in Fig. 21 The Fourier-transform confirms the 
dominant frequencies just stated, and reveals several 
more frequency components ranging up to 4kHz. The 
high spin evolutions S. 4z and S. 4x exhibit identical fre- 
qencies for all hve discernible peaks. In contrast, the 
low-spin evolution S-. 05z shows different frequencies. 

We interpret these features as a collection of excited 
quasi-normal modes in each neutron star. The modes are 
excited because the initial data is not precisely in equi¬ 
librium. For the two high-spin cases the neutron stars 
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FIG. 21. The Fourier transforms of the central density in all 
three of our runs. Labelled are the peak frequencies for the 
quasi-radial F mode and the I = 2, mode. 



FIG. 22. Fourier transforms of the central density pc{t), two 
modes of the magnitude of gravitational wave strain (|h 2 , 2 | 
and |h 2 ,o|), and x for the S.4z run. All quantities show 
excess power at 1.14kHz and 1.4kHz, corresponding to the 
frequencies of excited neutron star quasi-normal modes. 


have similar spin, and therefore the same quasi-normal 
modes, whereas in the low-spin model, the quasi-normal 
mode frequencies differ due to the different magnitude of 
the spin. 

To strengthen our interpretation, we consider the series 
of rotating, relativistic, T = 2 polytropes computed by 
Dimmelmeier et al [79] . 

Ref. [79] ’s model “AU3” has a central density of 
1.074 X 10“^Mq^ and its rotation is quantified through 
the ratio of polar to equatorial radius, r^/re = 0.780. 
Meanwhile, our high-spin runs have a central density 
of 1.02 X 10“^Mq^ (measured as time-average of the 
data shown in Fig. 201 and from our initial data, we 
find Tp/re ~ 0.8. Given the similarity in these values, 
we expect Ref. (Til’s “AU3” to approximate our high- 
spin stars S.4x, S.4z. Ref. m reports a frequency of 
fp = 1.283kHz for the spherically symmetric (£ = 0) 
F-mode, and a frequency f 2 f = 1.537kHz for the ax- 
isymmetric i = 2 mode ^/. These frequencies compare 
favorably with the two dominant frequencies in Fig. |21[ 
1.14kHz and 1.42kHz. 

Presumably, the small differences in these frequencies 
can be accounted for by the slight differences in stellar 
mass, radius, and rotation. Moreover, tidal interactions 
and orbital motion could factor in, as well. In our fig¬ 
ure]^ we also see several other peaks at higher frequen¬ 
cies, which are reminiscent of the overtones and mode 
couplings in figure 10 of m- If we identify our peak at 
fni = 4.03kHz with the Hi mode, then (in analogy to 
[fg Fig. 10), fni-fp = (4.03-1.14)kHz = 2.89kHz, and 
2fp = 2.28kHz, two frequencies that are indeed present 
in our simulations. Although we find clear indications 
of axisymmetric i = 2-modes, we note that their power 
is smaller by two orders of magnitude, compared to the 
spherically symmetric, dominant F mode. 

Turning to the low-spin run S.05z, we note that if, 


to first order, these frequencies scale like / ~ y/p (on 
dimensional grounds), then we expect to see F = 1.22kHz 
and = 1.49kHz. This is very close to what is seen. 

The density oscillations discussed in this section are re¬ 
flected in analogous oscillations in various other diagnos¬ 
tic quantities, for instance, the orbital frequency, Fig.fT^ 
and the quasi-local spin as shown in Fig. |14[ The dom¬ 
inant frequencies 1.14kHz and 1.42kHz can be robustly 
identified throughout our data analysis. In figure we 
plot the Fourier transform of the density, the (2,0) and 
(2, 2) gravitational wave strains, the orbital angular ve¬ 
locity time derivative dQ/dt and the measured spin x 
for the S. 4z run. All show peaks in power at these two 
frequencies, F ^ 1.14kHz and ~ 1.4kHz. Gold et. 
al. |80| . in the simulation of close encounters in eccen¬ 
tric, irrotational, NSNS binaries, find an excited f-mode 
frequency of 1.586 kHz. 

We believe that the stellar modes are excited because 
the initial data are not in perfect equilibrium. We expect 
the quasi-equilibrium approximations that enter the ini¬ 
tial data formalism to become less valid at higher spins, 
consistent with our observation that the high spin mod¬ 
els exhibit stronger oscillations. This interpretation is 
strengthened by additional simulations of neutron stars 
at larger separation. Increasing the intial separation by a 
factor 1.5, while keeping the same rotation parameter lo 
as in the S. 4z-case, we find quasi-normal oscillations of 
similar amplitude than in S.4z. If the oscillations were 
caused by the neglect of tidal deformation, we would ex¬ 
pect the amplitude to drop with the 3rd power of sepa¬ 
ration, inconsistent with our results. 


Finally, we point out that the radial rotation profile, cf. 
Eq. (48) influences the amplitude of the induced quasi¬ 
normal oscillations. If the initial data is constructed with 
the rotation profile Eq. (49), instead of equation 48 then 
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the amplitude of the density oscillations for high spin 
doubles. This further supports our conjecture that the 
origin of this mode comes from non-equilibrium initial 
data. 


V. DISCUSSION 

In this paper we implement Tichy’s method |18j to 
construct binary neutron star initial data with arbitrary 
rotation rates. We demonstrate that our implementation 
is exponentially convergent, as expected for the employed 
spectral methods. 

We measure the spin of the resulting neutron stars us¬ 
ing the quasi-local angular momentum formalism |40L 1431 
El [HI]. The resulting angular momentum is found to be 
nearly independent on the precise choice of extraction 
sphere, cf. Fig. and provides a means to define the 
quasi-local angular momentum of each neutron star to 
about 1%, both in the initial data and during the evolu¬ 
tion, cf. Fig.[^ We are able to construct binary neutron 
star initial data with dimensionless angular momentum 
of each star as large as y = S/M"^ ~ 0.43, both for the 
case of aligned spins, and also for a processing binary 
where the initial neutron star spins are tangential to the 
initial orbital plane. 

For irrotational BNS initial data sets, we find a quasi¬ 
local angular momentum of x ~ 2 x 10“^, cf. Fig. 
This spin is small enough that present waveform model¬ 
ing studies for BNS (e.g. [55H55] 1 are not yet limited by 
residual spin. 

When evolving the initial data sets, the dimensionless 
spin measured in the initial data drops by about 0.004, 
and then remains constant through the 10 inspiral orbits 
for which we evolved the neutron star binaries. During 
these evolutions, we also demonstrated iterative eccen¬ 
tricity removal: By analyzing the orbital frequency n{t) 
during the first few orbits, we can correct the initial data 
parameters Dg and dg, and thus decrease the orbital ec¬ 
centricity from e « 0.01 to e < 0.001. 

For the precessing simulation S.4x, we find precession 
of the neutron star spin directions. The numerically es¬ 
tablished precession of the spin axes and of the orbital 
angular momentum agrees well with post-Newtonian pre¬ 
dictions. 

The rotating neutron stars constructed here exhibit 
clear signals of exciting quasi-normal modes. We are 
able to identify multiple modes in the Fourier spectrum 
of the central density. The amplitude of the excited 
quasi-normal modes increases steeply with rotation 
rate of the neutron stars. For S-. 05z (spin magnitude 
X = 0.045) the density oscillations have peak-to-peak 
amplitude of 0.6%, raising to 20% for the two runs with 
high spins (S.4x and S.4z). 
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